#!/bin/sh

###############################################################################
##
##  HIGH VERBOSITY EXAMPLE
##
###############################################################################

# run from directory where this script is
cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
EXAMPLE_DIR=`pwd`

# check whether echo has the -e option
if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi

$ECHO
$ECHO "$EXAMPLE_DIR : starting"
$ECHO
$ECHO "This example shows how to calculate electron-phonon interaction"
$ECHO "coefficients for fcc Al."

# set the needed environment variables
. ../../../environment_variables

# required executables and pseudopotentials
BIN_LIST="pw.x ph.x q2r.x matdyn.x lambda.x"
PSEUDO_LIST="Al.pz-vbc.UPF"

$ECHO
$ECHO "  executables directory: $BIN_DIR"
$ECHO "  pseudo directory:      $PSEUDO_DIR"
$ECHO "  temporary directory:   $TMP_DIR"
$ECHO "  checking that needed directories and files exist...\c"

# check for directories
for DIR in "$BIN_DIR" "$PSEUDO_DIR" ; do
    if test ! -d $DIR ; then
        $ECHO
        $ECHO "ERROR: $DIR not existent or not a directory"
        $ECHO "Aborting"
        exit 1
    fi
done
for DIR in "$TMP_DIR" "$EXAMPLE_DIR/results" ; do
    if test ! -d $DIR ; then
        mkdir $DIR
    fi
done
cd $EXAMPLE_DIR/results

# check for executables
for FILE in $BIN_LIST ; do
    if test ! -x $BIN_DIR/$FILE ; then
        $ECHO
        $ECHO "ERROR: $BIN_DIR/$FILE not existent or not executable"
        $ECHO "Aborting"
        exit 1
    fi
done

# check for pseudopotentials
for FILE in $PSEUDO_LIST ; do
    if test ! -r $PSEUDO_DIR/$FILE ; then
       $ECHO
       $ECHO "Downloading $FILE to $PSEUDO_DIR...\c"
            $WGET $PSEUDO_DIR/$FILE \
                http://www.quantum-espresso.org/pseudo/1.3/UPF/$FILE 2> /dev/null
    fi
    if test $? != 0; then
        $ECHO
        $ECHO "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
        $ECHO "Aborting"
        exit 1
    fi
done
$ECHO " done"

# how to run executables
PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
PH_COMMAND="$PARA_PREFIX $BIN_DIR/ph.x $PARA_POSTFIX"
Q2R_COMMAND="$PARA_PREFIX $BIN_DIR/q2r.x $PARA_POSTFIX"
MATDYN_COMMAND="$PARA_PREFIX $BIN_DIR/matdyn.x $PARA_POSTFIX"
LAMBDA_COMMAND="$BIN_DIR/lambda.x "
$ECHO
$ECHO "  running pw.x as: $PW_COMMAND"
$ECHO "  running ph.x as: $PH_COMMAND"
$ECHO "  running q2r.x as: $Q2R_COMMAND"
$ECHO "  running matdyn.x as: $MATDYN_COMMAND"
$ECHO "  running lambda.x as: $LAMBDA_COMMAND"
$ECHO

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/*
$ECHO "  done"

#
# SCF at dense k-mesh, good enough for electronic DOS
#
cat > al.scf.fit.xml << EOF
<?xml version="1.0" encoding="UTF-8"?>


<input calculation="scf" prefix="al">

	<cell type="qecell">
		<qecell ibrav="2" alat="7.5">
			<real rank="1" n1="5">
				0.0 0.0 0.0 0.0 0.0
			</real>
		</qecell>
	</cell>

	<atomic_species ntyp="1">
		<specie name="Al">
			<property name="mass">
				<real>26.98</real>
			</property>
			<property name="pseudofile">
				<string>Al.pz-vbc.UPF</string>
			</property>
		</specie>
	</atomic_species>

	<atomic_list units="alat" nat="1" >
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.00 0.00 0.00
				</real>
			</position>
		</atom>						
	</atomic_list>	
	
	
	<field name="InputOutput">

		<parameter name="restart_mode">
			<string>
				from_scratch
			</string>
		</parameter>

		<parameter name="pseudo_dir">
			<string>
				$PSEUDO_DIR/
			</string>
		</parameter>
		
		<parameter name="outdir">
			<string>
				$TMP_DIR/
			</string>
		</parameter>

		<parameter name="la2F">
			<logical>
				true
			</logical>
		</parameter>
		
	</field>
	
	<field name="Numerics">

		<parameter name="ecutwfc">
			<real>
				15.0
			</real>
		</parameter>
		
		<parameter name="mixing_beta">
			<real>
				0.7
			</real>
		</parameter>
		
		<parameter name="conv_thr">
			<real>
				1.0d-8
			</real>
		</parameter>
	
	</field>
	
	<field name="Options">

		<parameter name="occupations">
			<string>
				smearing
			</string>
		</parameter>
		
		<parameter name="smearing">
			<string>
				methfessel-paxton
			</string>
		</parameter>
		
		<parameter name="degauss">
			<real>
				0.05
			</real>
		</parameter>		
	
	</field>	
	
	<k_points type="automatic">
		<mesh>
			<integer rank="1" n1="6">
				16 16 16  0 0 0
			</integer>
		</mesh> 
	</k_points>
</input>
EOF
$ECHO "  running the scf calculation with dense k-point grid...\c"
$PW_COMMAND  < al.scf.fit.xml > al.scf.fit.out
check_failure $?
$ECHO "  done"
#
#  SCF at k-mesh good enough for phonons
#
cat > al.scf.xml << EOF
<?xml version="1.0" encoding="UTF-8"?>


<input calculation="scf" prefix="al">

	<cell type="qecell">
		<qecell ibrav="2" alat="7.5">
			<real rank="1" n1="5">
				0.0 0.0 0.0 0.0 0.0
			</real>
		</qecell>
	</cell>

	<atomic_species ntyp="1">
		<specie name="Al">
			<property name="mass">
				<real>26.98</real>
			</property>
			<property name="pseudofile">
				<string>Al.pz-vbc.UPF</string>
			</property>
		</specie>
	</atomic_species>

	<atomic_list units="alat" nat="1" >
		<atom name="Al">
			<position>
				<real rank="1" n1="3">
					0.00 0.00 0.00
				</real>
			</position>
		</atom>						
	</atomic_list>		
	
	
	<field name="InputOutput">

		<parameter name="restart_mode">
			<string>
				from_scratch
			</string>
		</parameter>

		<parameter name="pseudo_dir">
			<string>
				$PSEUDO_DIR/
			</string>
		</parameter>
		
		<parameter name="outdir">
			<string>
				$TMP_DIR/
			</string>
		</parameter>
		
	</field>
	
	<field name="Numerics">

		<parameter name="ecutwfc">
			<real>
				15.0
			</real>
		</parameter>
		
		<parameter name="mixing_beta">
			<real>
				0.7
			</real>
		</parameter>
		
		<parameter name="conv_thr">
			<real>
				1.0d-8
			</real>
		</parameter>
	
	</field>
	
	<field name="Options">

		<parameter name="occupations">
			<string>
				smearing
			</string>
		</parameter>
		
		<parameter name="smearing">
			<string>
				methfessel-paxton
			</string>
		</parameter>
		
		<parameter name="degauss">
			<real>
				0.05
			</real>
		</parameter>		
	
	</field>	
	
	<k_points type="automatic">
		<mesh>
			<integer rank="1" n1="6">
				8 8 8  0 0 0
			</integer>
		</mesh> 
	</k_points>
</input>
EOF
$ECHO "  running the scf calculation...\c"
$PW_COMMAND < al.scf.xml > al.scf.out
check_failure $?
$ECHO "  done"
#
cat > al.elph.in << EOF
Electron-phonon coefficients for Al
 &inputph
  tr2_ph=1.0d-10,
  prefix='al',
  fildvscf='aldv',
  amass(1)=26.98,
  outdir='$TMP_DIR/',
  fildyn='al.dyn',
  electron_phonon='interpolated',
  trans=.true.,
  ldisp=.true.
  nq1=4, nq2=4, nq3=4
 /
EOF
$ECHO "  running the el-ph calculation...\c"
$PH_COMMAND < al.elph.in > al.elph.out
check_failure $?
$ECHO "  done"
#
#   q2r and matdyn
#
cat > q2r.in << EOF
 &input
  zasr='simple',  fildyn='al.dyn', flfrc='Al444.fc', la2F=.true.
 /
EOF
$ECHO "  running q2r...\c"
$Q2R_COMMAND < q2r.in > q2r.out
check_failure $?
$ECHO "  done"
#
#
#
cat > matdyn.in.freq << EOF
 &input
    asr='simple',  amass(1)=26.98,
    flfrc='Al444.fc', flfrq='Al444.freq', la2F=.true., dos=.false.
 /
  19
  0.000 0.0 0.0     0.0
  0.125 0.0 0.0     0.0
  0.250 0.0 0.0     0.0
  0.375 0.0 0.0     0.0
  0.500 0.0 0.0     0.0
  0.750 0.0 0.0     0.0
  1.000 0.0 0.0     0.0
  0.825 0.125 0.125 0.0
  0.750 0.250 0.250 0.0
  0.625 0.375 0.375 0.0
  0.500 0.500 0.500 0.0
  0.325 0.325 0.325 0.0
  0.250 0.250 0.250 0.0
  0.125 0.125 0.125 0.0
  0.000 0.000 0.000 0.0
  0.125 0.125 0.000 0.0
  0.250 0.250 0.000 0.0
  0.325 0.325 0.000 0.0
  0.500 0.500 0.000 0.0
EOF
$ECHO "  running matdyn for frequency calculation...\c"
$MATDYN_COMMAND < matdyn.in.freq > matdyn.out.freq
check_failure $?
$ECHO "  done"
#
#
#
cat > matdyn.in.dos << EOF
 &input
    asr='simple',  amass(1)=26.98,
    flfrc='Al444.fc', flfrq='Al444.freq', la2F=.true., dos=.true.
    fldos='phonon.dos', nk1=10, nk2=10, nk3=10, ndos=50
 /
EOF
$ECHO "  running matdyn for a2F(omega) calculation...\c"
$MATDYN_COMMAND < matdyn.in.dos > matdyn.out.dos
check_failure $?
$ECHO "  done"
#
# 
#
cat > lambda.in << EOF
10  0.12  1    ! emax (something more than highest phonon mode in THz), degauss, smearing method 
    8          ! Number of q-points for which EPC is calculated, 
    0.0000000  0.0000000  0.0000000   1.00  ! the first q-point, use kpoints.x program to calculate 
    -0.2500000 -0.2500000  0.2500000   8.00  ! q-points and their weight
    -0.5000000 -0.5000000  0.5000000   4.00  !
     0.0000000  0.0000000  0.5000000   6.00  ! 4th q-point, qx,qy,qz
    -0.2500000 -0.2500000  0.7500000  24.00  !
    -0.5000000 -0.5000000  1.0000000  12.00  !
     0.0000000  0.0000000  1.0000000   3.00  !
    -0.5000000  0.0000000  1.0000000   6.00  ! the last q-point 
elph. 0.000000. 0.000000. 0.000000 ! elph output file names,
elph.-0.250000. 0.250000.-0.250000 ! in the same order as the q-points before  
elph. 0.500000.-0.500000. 0.500000
elph. 0.000000. 0.500000. 0.000000
elph. 0.750000.-0.250000. 0.750000
elph. 0.500000. 0.000000. 0.500000
elph. 0.000000.-1.000000. 0.000000
elph.-0.500000.-1.000000. 0.000000
0.10                     ! \mu the Coloumb coefficient in the modified
                         ! Allen-Dynes formula for T_c (via \omega_log)
EOF
$ECHO "  running lambda.x for lambda calculation...\c"
$LAMBDA_COMMAND < lambda.in > lambda.out
check_failure $?
$ECHO "  done"




